!----------------------------------------------------------------------!
! Optional advection scheme for CCATT-BRAMS/BRAMS models version 4.2+  ! 
! Based on Walcek, 2000 (JGR) and Walcek and Aleksic, 1998 (ATENV).    ! 
! The scheme is highly conservative, monotonic and keeps mass mixing   ! 
! ratio positive definite. 					       ! 
! Implemented by Saulo Freitas (saulo.freitas@cptec.inpe.br) @ Jun/2009! 
!----------------------------------------------------------------------!

MODULE monotonic_adv
  
  USE node_mod, only:        &
          	      ibcon, &  !INTENT(IN)
          	      mynum, &  !INTENT(IN)
          	      mi0,   &  !INTENT(IN)
          	      mj0,   &  !INTENT(IN)
          	      mmyp,  &  !INTENT(IN)
          	      mmxp,  &  !INTENT(IN)
          	      mmzp   	!INTENT(IN)
  USE mem_grid, ONLY:        &
  	    	     dtlt,   & !INTENT(IN)
  	    	     ngrids, & !INTENT(IN)
  	    	     ngrid,  & !INTENT(IN)
	    	     dzt,    & !INTENT(IN)
	    	     dztn,   & !INTENT(IN)
            	     grid_g, & !INTENT(IN)	 
	    	     grid_g, & !INTENT(IN) 
	    	     naddsc, & !INTENT(IN) 
	             hw4       !INTENT(IN)	   
  
  USE mem_basic, ONLY: basic_g  !INTENT(IN) 
  
  USE micphys    ,  ONLY: level !INTENT(IN) 
  
  USE rconstants ,  ONLY: cp,p00,cv,rgas,cpi   !INTENT(IN)

  use mem_aer1, only : aerosol,num_scalar_aer_1st   !INTENT(IN)
  use mem_chem1, only : nspecies_transported        !INTENT(IN)

  use mem_scratch, only  : scratch  ! only scr1, inout
  
  use var_tables, only : scalar_tab & ! (var_p = IN, var_t = INOUT) 
                        ,num_scalar   ! (IN) 

  IMPLICIT NONE

  INTEGER , PARAMETER :: ON=1,OFF=0

  INTEGER , PARAMETER :: use_true_density  = 0 ! 0= OFF, 1=ON
  
  !- for theoretical experiments
  INTEGER , PARAMETER :: theor_wind = 0        ! 0= OFF, 1=ON

  INTEGER           :: mnt_adv_initialized=0
  
  real, parameter :: c1 = cv/rgas, c2 = p00/rgas !c2 = p00*(cpi**c1)/rgas

  TYPE advmnt_vars 
     REAL,POINTER,DIMENSION  (:,:,:)  :: u3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: v3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: w3d 

     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3du 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dv 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dw 

     REAL,POINTER,DIMENSION  (:,:,:)  :: den0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den1_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den2_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den3_3d 

     REAL,POINTER,DIMENSION  (:,:)    :: dxtW 
     REAL,POINTER,DIMENSION  (:,:)    :: dytW 
     REAL,POINTER,DIMENSION  (:)      :: dztW 

  END TYPE advmnt_vars

  TYPE(advmnt_vars), ALLOCATABLE,DIMENSION(:)   :: advmnt_g

  PUBLIC :: advmnt_driver ! Subroutine

CONTAINS
  !----------------------------------------------------

  SUBROUTINE advmnt_driver(varn,m1,m2,m3,ia,iz,ja,jz,izu,jzv)
    IMPLICIT NONE
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    INTEGER , INTENT(IN) :: izu
    INTEGER , INTENT(IN) :: jzv
    character(len=*),intent(IN) :: varn

    !--- local vars
    integer n,ng,mxyzp
    real, pointer :: scalarp, scalart
    
    !-scratch arrays initialization  and air density for using
    !- basic state densities
    IF(mnt_adv_initialized == OFF) THEN
       CALL initialize_advmnt(ngrids,mmzp,mmxp,mmyp)
       
       do ng=1,ngrids
         CALL initialize_grid_spacings(ng,mmzp(ng),mmxp(ng),mmyp(ng) &
                      ,grid_g(ng)%dxt	    &
                      ,grid_g(ng)%dyt	    &
 		      ,grid_g(ng)%fmapt     &
		      ,grid_g(ng)%rtgt      &
!
    		      ,advmnt_g(ng)%dxtW &
    		      ,advmnt_g(ng)%dytW &
		      ,advmnt_g(ng)%dztW )
       enddo
       
       if(use_true_density == OFF) then 
          do ng=1,ngrids
	    CALL initialize_densities(mmzp(ng),mmxp(ng),mmyp(ng) &
                      , basic_g(ng)%dn0     &
		      , basic_g(ng)%dn0u    &
                      , basic_g(ng)%dn0v    &
                      ,advmnt_g(ng)%dd0_3d  &
    		      ,advmnt_g(ng)%dd0_3du &
    		      ,advmnt_g(ng)%dd0_3dv &
		      ,advmnt_g(ng)%dd0_3dw )
	   enddo
       endif
       
       mnt_adv_initialized= ON    
    ENDIF
    mxyzp=m1*m2*m3
    
    !- Advect  U, V, and W
    if (varn .eq. 'V' .or. varn .eq. 'ALL') then
      stop 'not using mnt to advect u,v,w'
    endif

    !- Advect  scalars
    
    !- get actual air densities, if using them instead of basic state fields
    if(use_true_density == ON) then
     call get_true_densities(m1,m2,m3,level &
                      , basic_g(ngrid)%rtp     &
		      , basic_g(ngrid)%rv      &
                      , basic_g(ngrid)%pp      &
		      , basic_g(ngrid)%pi0     &
                      , basic_g(ngrid)%theta   &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )

    endif
 
 
    !- prepare wind velocities including map factors
    call prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,basic_g(ngrid)%uc  &
		      ,basic_g(ngrid)%up  &
                      ,basic_g(ngrid)%vc  &
                      ,basic_g(ngrid)%vp  &
		      ,basic_g(ngrid)%wc  &
                      ,basic_g(ngrid)%wp  &
!        
 		      ,grid_g(ngrid)%fmapui &
		      ,grid_g(ngrid)%fmapvi &
		      ,grid_g(ngrid)%rtgt   &
                      ,grid_g(ngrid)%rtgu   &
		      ,grid_g(ngrid)%rtgv   &
		      ,grid_g(ngrid)%f13t   &
                      ,grid_g(ngrid)%f23t   &
!
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  )
                                             
 
    if(theor_wind == on) then
 
        call prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    endif
    !- prepare Walcek's air densities
    
    call get_Walceks_densities(dtlt,m1,m2,m3 &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw &
                      ,advmnt_g(ngrid)%den0_3d &
    		      ,advmnt_g(ngrid)%den1_3d &
    		      ,advmnt_g(ngrid)%den2_3d &
		      ,advmnt_g(ngrid)%den3_3d &
    		      ,advmnt_g(ngrid)%dxtW &
    		      ,advmnt_g(ngrid)%dytW &
		      ,advmnt_g(ngrid)%dztW &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    )

     !- ready to do advection, loop over all scalars
     do n=1,num_scalar(ngrid)     
     if (scalar_tab(n,ngrid)%name .ne. 'COP') cycle
     
     print*,' doing=',n,scalar_tab(n,ngrid)%name     

!srf - somente para gases e aerossois
!     do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
!srf - somente para gases e aerossois



 !srf-falta incluir controle para shaved_eta e aerossois com sedimentacao


      scalarp => scalar_tab(n,ngrid)%var_p
      scalart => scalar_tab(n,ngrid)%var_t
      call atob(mxyzp,scalarp,scratch%scr1(1))

      call advect_mnt(m1,m2,m3,ia,iz,ja,jz,dtlt &
                     ,scalarp                 & ! initial mass mixing ratio (mmxr)
		     ,scratch%scr1(1)         & ! final mmxr after advection
                     ,advmnt_g(ngrid)%u3d     & ! 3D U wind
    		     ,advmnt_g(ngrid)%v3d     & ! 3D V wind
    		     ,advmnt_g(ngrid)%w3d     & ! 3D SIGMAZ wind
                     ,advmnt_g(ngrid)%dd0_3d  & ! density 0
                     ,advmnt_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
    		     ,advmnt_g(ngrid)%den1_3d & ! density D1
    		     ,advmnt_g(ngrid)%den2_3d & ! density D2
		     ,advmnt_g(ngrid)%den3_3d & ! density D3
    		     ,advmnt_g(ngrid)%dxtW    &
    		     ,advmnt_g(ngrid)%dytW    &
		     ,advmnt_g(ngrid)%dztW    &
                     ,grid_g(ngrid)%dxt       & ! true dx
                     ,grid_g(ngrid)%dyt       & ! true dy
                     ,dzt(1:m1) 	      & ! true dz
		     ,grid_g(ngrid)%rtgt      )
     
      call advtndc(m1,m2,m3,ia,iz,ja,jz        &
                  ,scalarp  ,scratch%scr1 (1)  &
                  ,scalart  ,dtlt,mynum        )

 
     enddo
!    stop 100
  END SUBROUTINE advmnt_driver
 !----------------------------------------------------
  
  SUBROUTINE initialize_advmnt(ngrids, mmzp,mmxp,mmyp)
   implicit none 
   INTEGER , INTENT(IN) :: ngrids
   INTEGER , INTENT(IN) :: mmxp(ngrids)
   INTEGER , INTENT(IN) :: mmyp(ngrids)
   INTEGER , INTENT(IN) :: mmzp(ngrids)
   
   INTEGER :: ng
   
   !maxz=maxval(mmzp(1:ngrids))
   !maxx=maxval(mmxp(1:ngrids))
   !maxy=maxval(mmyp(1:ngrids))
   
   
   ALLOCATE (advmnt_g(ngrids))

   do ng=1,ngrids
      ALLOCATE(advmnt_g(ng)%u3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%u3d=0.
      ALLOCATE(advmnt_g(ng)%v3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%v3d=0.
      ALLOCATE(advmnt_g(ng)%w3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%w3d=0.
      
      ALLOCATE(advmnt_g(ng)%dd0_3d (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3d =0.
      ALLOCATE(advmnt_g(ng)%dd0_3du(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3du=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dv(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dv=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dw(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmnt_g(ng)%den0_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den0_3d=0.
      ALLOCATE(advmnt_g(ng)%den1_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den1_3d=0.
      ALLOCATE(advmnt_g(ng)%den2_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den2_3d=0.
      ALLOCATE(advmnt_g(ng)%den3_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den3_3d=0.

      ALLOCATE(advmnt_g(ng)%dxtW(	  mmxp(ng),mmyp(ng))); advmnt_g(ng)%dxtW=0.
      ALLOCATE(advmnt_g(ng)%dytW(	  mmxp(ng),mmyp(ng))); advmnt_g(ng)%dytW=0.
      ALLOCATE(advmnt_g(ng)%dztW(mmzp(ng)))                  ; advmnt_g(ng)%dztW=0.
   enddo

   



  END SUBROUTINE initialize_advmnt
 !----------------------------------------------------
 
  SUBROUTINE initialize_densities(m1,m2,m3,dn0,dn0u,dn0v &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::dn0,dn0u,dn0v
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::dd0_3d,dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   
   dd0_3d (:,:,:)=  dn0 (:,:,:) 
   dd0_3du(:,:,:)=  dn0u(:,:,:)
   dd0_3dv(:,:,:)=  dn0v(:,:,:)
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dn0(k,i,j) +dn0(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo
   

  END SUBROUTINE initialize_densities
 !----------------------------------------------------
 
  SUBROUTINE initialize_grid_spacings(ng,m1,m2,m3,dxt,dyt,fmapt,rtgt &
			    ,dxtW,dytW,dztW )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: ng,m1,m2,m3
   REAL,DIMENSION(m2,m3),intent(IN) :: dxt,dyt,fmapt,rtgt
   
   REAL,DIMENSION(m2,m3),intent(OUT):: dxtW,dytW
   REAL,DIMENSION(m1),intent(OUT):: dztW
   ! local var
   integer i,j,k
   real rtgti
   
   do j = 1,m3
      do i = 1,m2
   	rtgti = 1. / rtgt(i,j)
        
	!- at init/rams_grid.f90:
        !     dxt(i,j)=fmapt(i,j)/(xmn(i,ngrid)-xmn(i-1,ngrid))  
        !     dyt(i,j)=fmapt(i,j)/(ymn(j,ngrid)-ymn(j-1,ngrid))

   	 dxtW(i,j) = 1./(dxt(i,j) * fmapt(i,j) * rtgti)
   	 dytW(i,j) = 1./(dyt(i,j) * fmapt(i,j) * rtgti)
   enddo;enddo
   do k = 1,m1
    !- at init/gridset.f90:	
    !  dztn(k,ifm) = 1. / (zmn(k,ifm) - zmn(k-1,ifm))
    ! Por que o Jacobiano nao depende de Z, o dztw depende somente 
    ! de z.
    !dztW(k,i,j) = 1./ ( dzt(k) * rtgti * fmapt(i,j)**2 )
     dztW(k)	 = 1./ ( dztn(k,ng) ) !
 
   enddo
  END SUBROUTINE initialize_grid_spacings
 !----------------------------------------------------
 
  SUBROUTINE get_true_densities(m1,m2,m3,level,rtp,rv,pp,pi0,theta &
			       ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,level
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::rtp,rv,pp,pi0,theta
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::dd0_3d,dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   real c3
   
   c3 = c2 * (cpi**c1)

   !- true air density at points "T"

   if( level == 0 ) then
     dd0_3d(:,:,:) = (c3/theta(:,:,:))*(pi0(:,:,:)+pp(:,:,:))**c1
   else
     dd0_3d(:,:,:) = (c3/theta(:,:,:))* (1. + rtp(:,:,:))/(1. + 1.61*rv(:,:,:))*(pi0(:,:,:)+pp(:,:,:))**c1
   endif 

   !- true air density at points "U", "V" and "W":
   
   call fill_dn0uv(m1,m2,m3,dd0_3d,dd0_3du,dd0_3dv)

   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dd0_3d(k,i,j) + dd0_3d(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo

 END SUBROUTINE get_true_densities

 !----------------------------------------------------

 SUBROUTINE prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,uc,up,vc,vp,wc,wp &
  		            ,fmapui &
		            ,fmapvi &
                            ,rtgt   &
                            ,rtgu   &
		            ,rtgv   &
		            ,f13t   &
                            ,f23t   &
			    ,u3d,v3d,w3d)
   
   implicit none 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL,DIMENSION(m1,m2,m3),intent(IN) :: uc,up,vc,vp,wc,wp
   REAL,DIMENSION(m2,m3)   ,intent(IN) :: rtgt,rtgu,rtgv,fmapui,fmapvi,f13t,f23t
   
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d
  
   !- local var
   real   dtlto2
   integer jm,jp,im,ip 
   integer i,j,k
   real :: c1,c2,rtgti

   dtlto2 = .5! * dtlt

  ! u3d, u3d, and w3d are input as the velocity components (averaged
  ! between past and current time levels) times dtlt.
   do j=1,m3
     do i = 1,m2
      do k = 1,m1
  
          w3d(k,i,j) = ( wc(k,i,j) + wp(k,i,j) )*dtlto2
          u3d(k,i,j) = ( uc(k,i,j) + up(k,i,j) )*dtlto2
          v3d(k,i,j) = ( vc(k,i,j) + vp(k,i,j) )*dtlto2
  
   enddo;enddo;enddo


  !return ! for pure cartesian coordinates

  ! here w3d is the cartesian vertical velocity 
  
  ! Add contribution to w3d from horiz winds crossing sloping sigma surfaces,
  ! and include 1/rtgt factor in w3d
  do j = 1,m3		 
     jm = max(1,j-1)
     jp = min(m3,j+1)
     do i = 1,m2
  	im = max(1,i-1)
  	ip = min(m2,i+1)
        rtgti = 1. / rtgt(i,j)
	
	do k = 1,m1-1
	    w3d(k,i,j) = ( (u3d(k,i,j) + u3d(k+1,i,j) + u3d(k,im,j) + u3d(k+1,im,j) ) * f13t(i,j)  &
	               +   (v3d(k,i,j) + v3d(k+1,i,j) + v3d(k,i,jm) + v3d(k+1,i,jm) ) * f23t(i,j)  &
		         ) * hw4(k)  &
	               + w3d(k,i,j) * rtgti

  	enddo
     enddo
  enddo
  ! here w3d is the sigma_z velocity 

  !- including map factors on U,V:
  do j = 1,m3
     do i = 1,m2
  	c1 = fmapui(i,j) * rtgu(i,j)
  	c2 = fmapvi(i,j) * rtgv(i,j)
  	do k = 1,m1-1
  	   u3d(k,i,j) = u3d(k,i,j) * c1 
  	   v3d(k,i,j) = v3d(k,i,j) * c2 
  	enddo
     enddo
  enddo

  END SUBROUTINE prepare_winds
 !----------------------------------------------------

  SUBROUTINE get_Walceks_densities(dt,m1,m2,m3,u3d,v3d,w3d &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw &
			    ,den0_3d,den1_3d,den2_3d,den3_3d &
			    ,dxtW,dytW,dztW,dxt,dyt)

   IMPLICIT NONE 
   !-in
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL    , INTENT(IN) :: dt
   REAL, DIMENSION(m1)      , INTENT(IN) :: dztW
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL, DIMENSION(m1,m2,m3), INTENT(IN) :: u3d,v3d,w3d     &
			                   ,dd0_3d ,dd0_3du &
					   ,dd0_3dv,dd0_3dw
   !-out
   REAL, DIMENSION(m1,m2,m3), INTENT(OUT):: den0_3d,den1_3d &
                                           ,den2_3d,den3_3d
   
   ! local var
   integer i,j,k
   
    DO  j=m3,2,-1
      DO  i=2,m2
        DO k = 2,m1
	
            den0_3d(k,i,j)=dd0_3d(k,i,j)
	    
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxt(i,j)*(d0         *u3d(k,i,j)- d0           *u3d(k,i-1,j ))
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt*dxt(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))
	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxtW(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))

!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dx*(d0*v3d(k,i,j)- d0*v3d(k,i,j-1 ))
!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt*dyt(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
            den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dytW(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
!
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dx         *(d0*w3d(k,i,j)-d0*w3d(  k-1,i,j))
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt*dzt(k)     *(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))
            den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dztW(k)*(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))

       END DO !i
      END DO !j
    END DO !k  

 

 END SUBROUTINE get_Walceks_densities
 !----------------------------------------------------
 SUBROUTINE advect_mnt(m1,m2,m3,ia,iz,ja,jz,dt,vc3din,vc3dout_Z& 
                      ,u3d,v3d,w3d,dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d &
		      ,dxtW,dytW,dztW,dxt,dyt,dzt,rtgt)
 
  IMPLICIT NONE
  INTEGER , INTENT(IN) :: m1
  INTEGER , INTENT(IN) :: m2
  INTEGER , INTENT(IN) :: m3
  INTEGER , INTENT(IN) :: ia
  INTEGER , INTENT(IN) :: iz
  INTEGER , INTENT(IN) :: ja
  INTEGER , INTENT(IN) :: jz
  REAL    , INTENT(IN) :: dt
  
  REAL,DIMENSION(m1)      , INTENT(IN)  :: dzt
  REAL,DIMENSION(m2,m3)   , INTENT(IN)  :: dxt,dyt,rtgt
  
  REAL,DIMENSION(m1),       INTENT(IN) :: dztW 
  REAL,DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

  REAL,DIMENSION(m1,m2,m3), INTENT(IN) :: vc3din    ! initial mixing ratio 
  REAL,DIMENSION(m1,m2,m3), INTENT(IN) :: u3d,v3d,w3d
  REAL,DIMENSION(m1,m2,m3), intent(IN) :: dd0_3d,den0_3d,den1_3d,den2_3d,den3_3d

  REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3dout_Z ! updated mixing ratio due advection

  !- local var
  REAL,DIMENSION(m1)               :: dxx  
  REAL,DIMENSION(m2,m3)            :: dxy  
  real masscon,initialmass,vol
  integer :: i,j,k,nrec,it

!---------------------------------------------------------------------  
!- for theoretical experiments
!
!  OPEN(19,file='walcek_mixingratio_3d_rot_wind.gra',         &
!       form='unformatted',access='direct',status='unknown',  &
!       recl=4*m2*m3)
!  nrec =0
!  do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) u3d(k,:,:)
!  enddo
!  do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) v3d(k,:,:)
!  enddo
!  CLOSE(19)
!
!  OPEN(19,file='walcek_mixingratio_3d_rot.gra',         &
!       form='unformatted',access='direct',status='unknown',  &
!       recl=4*m2*m3)
!  nrec =0
!  
!  DO  it=1,3500
!
!   do k=1,m1
!    nrec=nrec+1
!    write (19,rec=nrec) vc3din(k,:,:)
!   enddo
! !-------------------------------
!  !- check mass conservation (initial mass)
!  initialmass=0.
!  do j=2,m3-1; do i=2,m2-1
!    do k=2,m1-1
!    !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
!    vol= 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j) ! true volume
!    initialmass=initialmass+vc3din(k,i,j)*dd0_3d(k,i,j)*vol
!
!  enddo;enddo;enddo
!---------------------------------------------------------------------
!---------------------------------------------------------------------
  
  ! Do I-advection first over all J-k lines
  !dxy(:,:)=1./dxt(:,:)
  !call Advec3d_X(m1,m2,m3,ia,iz,ja,jz,vc3din,u3d,den0_3d, den1_3d,dt,dxy,dd0_3d,vc3din)
!orig   call Advec3d_X(m1,m2,m3,ia,iz,ja,jz,vc3din,u3d,den0_3d, den1_3d,dt,dxtW,dd0_3d,vc3din)

  !- do J-advection next over all I-k lines
  !dxy(:,:)=1./dyt(:,:)
  !call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3din,v3d,den1_3d, den2_3d,dt,dxy,dd0_3d,vc3din)
!orig  call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3din,v3d,den1_3d, den2_3d,dt,dytW,dd0_3d,vc3din)
  call Advec3d_Y(m1,m2,m3,ia,iz,ja,jz,vc3din,v3d,den1_3d, den2_3d,dt,dytW,dd0_3d,vc3dout_Z)

  !- do k-advection next over all I-J lines
  ! vc3dout_Z = razao de mistura final (apos a adveccao nas 3 dimensoes)
  !dxx(1:m1)=1./dzt(1:m1)!tmp (usa dxx e nao dxy)
  !call Advec3d_Z(m1,m2,m3,ia,iz,ja,jz,vc3din,w3d,den2_3d, den3_3d,dt,dxx,dd0_3d,vc3dout_Z)
!orig  call Advec3d_Z(m1,m2,m3,ia,iz,ja,jz,vc3din,w3d,den2_3d, den3_3d,dt,dztW,dd0_3d,vc3dout_Z)


 
  !-------------------------------
  !- check mass conservation (after advection)
  !masscon=0.
  !do j=2,m3-1; do i=2,m2-1
  !  do k=2,m1-1
  !  !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
  !  vol= 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j) ! true volume
  !  masscon=masscon+vc3dout_Z(k,i,j)*dd0_3d(k,i,j)*vol
  !
  !enddo;enddo;enddo
  !print*,'masscon=',initialmass,masscon,100.*(masscon-initialmass)/initialmass
  !print*,'MAXVALS=', maxval(vc3din)*28.96/28, maxval(vc3dout_Z)*28.96/28
  !call flush(6)
  !-------------------------------

!---------------------------------------------------------------------
!- for theoretical experiments
!---------------------------------------------------------------------
!  vc3din=vc3dout_Z
!
! ENDDO
! close (19)
! stop 333
!---------------------------------------------------------------------
!---------------------------------------------------------------------




 END SUBROUTINE advect_mnt
 !----------------------------------------------------
 
SUBROUTINE Advec3d_X(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m2)                 :: flux
  REAL,DIMENSION(m2)                   :: vcmax
  REAL,DIMENSION(m2)                   :: vcmin
  LOGICAL,DIMENSION(m2)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

 idime = m2 ! x dir
 !- big loop y-z
 DO k=2,m1 ; DO j=ja,jz
 
  
  imxmn=.false.
 ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1 ! ia,iz
  
     imxmn(i)=q0(k,i,j)>=max(q0(k,i-1,j),q0(k,i+1,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i-1,j),q0(k,i+1,j))       !       extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k,i,j  )< zr0) ck1= q0(k,i+1,j)
     if(u(k,i-1,j)>=zr0) ck2= q0(k,i-1,j)
     vcmax(i)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(i)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(k,1,j)>=zr0) flux(1)= q0(k,1,j)*u(k,1,j)*dt*dd0(k,1,j)
  
  DO i=2,idime-1 ! ia,iz
  
  !print*,'part 1'
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k,i-1,j)<zr0) THEN
        flux(i)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream
        x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
        x1n= (1.-x1)*(q0(k,i+1,j)-q0(k,i-1,j))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(i-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(i+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i+1,j))  ), MAX(q0(k,i,j),q0(k,i+1,j)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(i),   MIN(   vcmax(i),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(i-1)/dxx(i,j))/den1(k,i,j) ))

	! for debug purposes only
        if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(I-1)/DXX(I,J))/DEN1(k,i,j)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(i)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(i-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  
  IF(u(k,idime-1,j)<zr0) flux(idime-1)=q0(k,idime,j)*u(k,idime-1,j)*dt*dd0(k,idime-1,j)
 
  DO i=idime-1,2,-1 !iz,ia,-1
 ! print*,'part 2'
     
       IF(u(k,i-1,j)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(i),   MIN(   vcmax(i),&
           (q0(k,i,j)*den0(k,i,j)-flux(i)/dxx(i,j) + flux(i-1)/dxx(i,j))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k,i-1,j))/dxx(i,j)     ! Courant number
           x1n= (1.-x1)*(q0(k,i-1,j)-q0(k,i+1,j))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(i+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(i-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i-1,j)) ), MAX(q0(k,i,j),q0(k,i-1,j)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(i),  MIN(   vcmax(i), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(i)/dxx(i,j)-x1*cf1*dd0(k,i-1,j))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(I)/DXX(I,J)-X1*CF1*DD0(k,I-1,j))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(i-1)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(i)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop y-z

END SUBROUTINE Advec3d_X

!---------------------------------------------------------------------------

SUBROUTINE Advec3d_Y(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(inOUT),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(inOUT),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(inOUT),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(inOUT),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(inOUT)                   :: dt
  REAL   ,INTENT(inOUT),DIMENSION(m2,m3)     :: dxx
  REAL   ,INTENT(inOUT),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m3)                 :: flux
  REAL,DIMENSION(m3)                   :: vcmax
  REAL,DIMENSION(m3)                   :: vcmin
  LOGICAL,DIMENSION(m3)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

!u=1.0
!dd0=1.
!den1=1.;den0=1.
!dxx=1000.
!dt = 1.

 idime = m3 ! y dir
 !- big loop y-z
 DO k=10,10 ; DO i=50,50
!  flux(1:m3)=0.
  
  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  j=2,idime-1 ! ja,jz
  
     imxmn(j)=q0(k,i,j)>=max(q0(k,i,j-1),q0(k,i,j+1)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i,j-1),q0(k,i,j+1))	    !	    extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k,i,j  )< zr0) ck1= q0(k,i,j+1)
     if(u(k,i,j-1)>=zr0) ck2= q0(k,i,j-1)
     vcmax(j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  
  go to 100
  
  IF(u(k,i,1)>=zr0) flux(1)= q0(k,i,1)*u(k,i,1)*dt*dd0(k,i,1)
  
  DO j=2,idime-1 ! ja,jz
 
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k,i,j-1)<zr0) THEN
        flux(j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream


        x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
        x1n= (1.-x1)*(q0(k,i,j+1)-q0(k,i,j-1))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
	!write(mynum+1,*) '1' ,j,flux(j),q0(k,i,j)!,x1,x1n;
	!call flush(mynum+1)
	
	
	
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(j-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(j+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j+1))  ), MAX(q0(k,i,j),q0(k,i,j+1)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(j),   MIN(   vcmax(j),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(j-1)/dxx(i,j))/den1(k,i,j) ))


!	write(mynum+1,*) '2', j,qn(k,i,j),cf1,cf,q0(k,i,j+1);call flush(mynum+1)
	write(mynum+1,*) '2', j,qn(k,i,j) ,q0(k,i,j-1),flux(j-1)
	if(j-1 .eq.1 .or. j-1.eq.65) write(mynum+1,*) 'X', j,qn(k,i,j-1);call flush(mynum+1)


        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(j-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  cycle!XXXXXXXXXXXXXX
  100 continue
  
  IF(u(k,i,idime-1)<zr0) flux(idime-1)=q0(k,i,idime)*u(k,i,idime-1)*dt*dd0(k,i,idime-1)
 
  DO j=idime-1,2,-1 !jz,ja,-1
     
       IF(u(k,i,j-1)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(j),   MIN(   vcmax(j),&
           (q0(k,i,j)*den0(k,i,j)-flux(j)/dxx(i,j) + flux(j-1)/dxx(i,j))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k,i,j-1))/dxx(i,j)     ! Courant number
           x1n= (1.-x1)*(q0(k,i,j-1)-q0(k,i,j+1))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(j+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(j-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j-1)) ), MAX(q0(k,i,j),q0(k,i,j-1)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(j),  MIN(   vcmax(j), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(j)/dxx(i,j)-x1*cf1*dd0(k,i,j-1))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(J)/DXX(I,J)-X1*CF1*DD0(k,I,j-1))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(j-1)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(j)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop x-z


END SUBROUTINE Advec3d_Y
!---------------------------------------------------------------------------
SUBROUTINE Advec3d_Z(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m1)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1)                 :: flux
  REAL,DIMENSION(m1)                   :: vcmax
  REAL,DIMENSION(m1)                   :: vcmin
  LOGICAL,DIMENSION(m1)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

 idime = m1 ! z dir
 !- big loop y-x
 DO j=ja,jz ; DO i=ia,iz
 
  
  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  k=2,idime-1 !
  
     imxmn(k)=q0(k,i,j)>=max(q0(k-1,i,j),q0(k+1,i,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k-1,i,j),q0(k+1,i,j))	    !	    extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k  ,i,j)< zr0) ck1= q0(k+1,i,j)
     if(u(k-1,i,j)>=zr0) ck2= q0(k-1,i,j)
     vcmax(k)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(k)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1,i,j)>=zr0) flux(1)= q0(1,i,j)*u(1,i,j)*dt*dd0(1,i,j)
  
  DO k=2,idime-1 ! ja,jz
  
  !print*,'part 1'
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k-1,i,j)<zr0) THEN
        flux(k)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream
        x1= dt*u(k,i,j)/dxx(k)               ! Courant number
        x1n= (1.-x1)*(q0(k+1,i,j)-q0(k-1,i,j))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(k-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(k+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k+1,i,j))  ), MAX(q0(k,i,j),q0(k+1,i,j)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(k),   MIN(   vcmax(k),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k-1)/dxx(k))/den1(k,i,j) ))

	! for debug purposes only
        if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k-1)/dxx(k))/DEN1(k,i,j)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(k)= dxx(k)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  
  IF(u(idime-1,i,j)<zr0) flux(idime-1)=q0(idime,i,j)*u(idime-1,k,j)*dt*dd0(idime-1,i,j)
 
  DO k=idime-1,2,-1 !jz,ja,-1
     
       IF(u(k-1,i,j)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k),   MIN(   vcmax(k),&
           (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k) + flux(k-1)/dxx(k))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k-1,i,j))/dxx(k)     ! Courant number
           x1n= (1.-x1)*(q0(k-1,i,j)-q0(k+1,i,j))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(k+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k-1,i,j)) ), MAX(q0(k,i,j),q0(k-1,i,j)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(k),  MIN(   vcmax(k), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k)-x1*cf1*dd0(k-1,i,j))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k)/dxx(k)-X1*CF1*DD0(k-1,I,j))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(k-1)=dxx(k)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop y-x


END SUBROUTINE Advec3d_Z
!-------------------------
SUBROUTINE Advec1d(idime,ix,q0,u,den0,den1,dt,dxx,dd0,qn,line)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: ix,line
  INTEGER,INTENT(IN)                   :: idime
  REAL   ,INTENT(in),DIMENSION(ix) :: q0
  REAL   ,INTENT(in),DIMENSION(ix)   :: u
  REAL   ,INTENT(in),DIMENSION(ix)     :: den0
  REAL   ,INTENT(in),DIMENSION(ix)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(ix)     :: dxx
  REAL   ,INTENT(in),DIMENSION(ix)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(ix)    :: qn

  REAL,DIMENSION(ix)                 :: flux
  REAL,DIMENSION(ix)                   :: vcmax
  REAL,DIMENSION(ix)                   :: vcmin
  LOGICAL,DIMENSION(ix)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1
     imxmn(i)=q0(i)>=max(q0(i-1),q0(i+1)) .or. & !=true if local
              q0(i)<=min(q0(i-1),q0(i+1))       !       extrema
     ck1= q0(i)
     ck2= q0(i)
     if(u(i  )< zr0) ck1= q0(i+1)
     if(u(i-1)>=zr0) ck2= q0(i-1)
     vcmax(i)= max( q0(i), ck1, ck2 )                      ! Eq-7
     vcmin(i)= min( q0(i), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1)>=zr0) flux(1)= q0(1)*u(1)*dt*dd0(1)
  DO i=2,idime-1
  
  !print*,'part 1'
     IF(u(i)<zr0) CYCLE
     IF(u(i-1)<zr0) THEN
        flux(i)= q0(i)*u(i)*dt*dd0(i)    !  outflow-only cell
     ELSE                              !      use upstream
        x1= dt*u(i)/dxx(i)               ! Courant number
        x1n= (1.-x1)*(q0(i+1)-q0(i-1))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(i) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(i-1)) cf= q0(i) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(i+1)) cf= q0(i) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(I) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(i),q0(i+1))  ), MAX(q0(i),q0(i+1)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(i)= MAX(  vcmin(i),   MIN(   vcmax(i),          &   !eq-3&8
        (q0(i)*den0(i)-x1*cf1*dd0(i)+flux(i-1)/dxx(i))/den1(i) ))

	! for debug purposes only
        if(debug==on)  QN(I)=(Q0(I)*DEN0(I)-X1*CF1*DD0(I)+FLUX(I-1)/DXX(I))/DEN1(I)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(i)= dxx(i)*(q0(i)*den0(i) - qn(i)*den1(i)) + flux(i-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  IF(u(idime-1)<zr0) flux(idime-1)=q0(idime)*u(idime-1)*dt*dd0(idime-1)
  DO i=idime-1,2,-1
 ! print*,'part 2'
       IF(u(i-1)>=zr0) THEN           ! Inflow-only cell
        IF(u(i)<zr0) qn(i)=  MAX(  vcmin(i),   MIN(   vcmax(i),&
           (q0(i)*den0(i)-flux(i)/dxx(i) + flux(i-1)/dxx(i))/den1(i) ))
        ELSE
           x1=  dt*ABS(u(i-1))/dxx(i)     ! Courant number
           x1n= (1.-x1)*(q0(i-1)-q0(i+1))/4.
           cf= q0(i) + x1n                                       !Eq-4b
           IF(imxmn(i+1)) cf= q0(i) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(i-1)) cf= q0(i) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(i),q0(i-1)) ), MAX(q0(i),q0(i-1)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(i)>=zr0) cf1= q0(i)     ! outflow-only cell upstream
           qn(i)= MAX(  vcmin(i),  MIN(   vcmax(i),           &   !Eq-3&8
               (q0(i)*den0(i)-flux(i)/dxx(i)-x1*cf1*dd0(i-1))/den1(i) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(I)=(Q0(I)*DEN0(I)-FLUX(I)/DXX(I)-X1*CF1*DD0(I-1))/DEN1(I)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(i-1)=dxx(i)*(qn(i)*den1(i) - q0(i)*den0(i)) + flux(i)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
  END SUBROUTINE Advec1d
!-------------------------------------------------------------------
 
  SUBROUTINE prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,u3d,v3d,w3d&
			    ,dxt,dyt &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )

   
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d     &
                                        ,dd0_3d ,dd0_3du &
					,dd0_3dv,dd0_3dw 
  
   !- local var
   REAL   :: dtlto2
   INTEGER :: i,j,k
   REAL  :: ai0s  =  25.0
   REAL  :: aj0s  =  50.0
   REAL  :: umx   =   10.0
   REAL,PARAMETER :: pii   =   3.141592653589793
   REAL    :: umax  =   0.0
   REAL    :: anrev,curnt,rx,xa,ilop,iwndty,nrec
   
   dtlto2 =  10.!*dtlt
    
     !WRITE(6,*) ' Wind fields?  0-rotating; or  1-divergent winds'
     iwndty = 0  ! 0-rotating
     !iwndty = 1  ! 1-divergent winds

     IF(iwndty==1) ai0s= 50.5
     ilop= ai0s-21.  ! needed for printouts
     ! Define wind fields (rotation or divergent) and initial mixing ratios
     !  Cone at (25,50) for rotating winds; Cone at (50,50) divergent winds
     DO k = 1,m1
   	DO  j=m3,1,-1
   	   DO  i=1,m2
   	   
   	       dd0_3d (k,i,j)=1.
	       dd0_3du(k,i,j)=1. 
               dd0_3dv(k,i,j)=1.
	       dd0_3dw(k,i,j)=1. 

   	       u3d(k,i,j)= -2.*umx*(REAL(j)-REAL(110)/2.-.5)/REAL(110)*dtlto2
   	       v3d(k,i,j)=  2.*umx*(REAL(i)-REAL(100)/2.-.5)/REAL(100)*dtlto2
 	       w3d(k,i,j)= 0.                                        *dtlto2
   	   
   	      IF(iwndty==1) THEN
   		 xa=pii/25.
   		 IF(J>0) u3d(k,i,j)=umx*SIN(xa*REAL(i))*SIN(xa*(REAL(j)))
   		 IF(I>0) v3d(k,i,j)=umx*COS(xa*(REAL(i)-.5))*COS(XA*(REAL(j)+.5))
   	      END IF
   	   
   	   
   	      umax= MAX(ABS(u3d(k,i,j)),ABS(v3d(k,i,j)),umax)
   	      rx= SQRT((REAL(i)-ai0s)**2.+(REAL(j)-aj0s)**2.)

   	   END DO !i
   	END DO !j
     END DO !k  

     !IF(iwndty==0) THEN
   	!WRITE(6,*) ' #iterations per revolution? must be>',pii*(m2-1)
   	!WRITE(6,*) '	Followed by the number of revolutions desired'
   	!srf- READ(5,*) NITPREV, ANREV
   	!nitprev=320; anrev=1
   	!dt= pii*REAL(m2-1)*dx/umax/REAL(nitprev)
   	!nstep= REAL(nitprev)*anrev + .9999

     !ELSE
   	!WRITE(6,*) ' Courant #?  eg. C=.6594 (to match Staniforth figs)'
   	!WRITE(6,*) '  Followed by # of iterations? e.g. 20,40,60,80,400'
   	! READ(5,*) CURNT, NSTEP
   	!curnt=0.6594 ;nstep=80
   	!dt= curnt*dx/umax
     !END IF

     !WRITE(6,*) ' Max wind speed, Courant number=',umax,umax*dt/dx
     !WRITE(6,*) nstep,' iterations required'
     !CALL FLUSH(6)


  END SUBROUTINE prepare_theor_winds
 !----------------------------------------------------------------------------  

  SUBROUTINE check_mass_conservation
  use mem_grid
  use var_tables
  use node_mod
  use mem_basic
  use mem_aer1, only : AEROSOL,num_scalar_aer_1st
  use mem_chem1, only : NSPECIES_TRANSPORTED
  
  implicit none
  real masscon, initalmass
  integer :: mxyzp,n
  mxyzp = mxp * myp * mzp
  initalmass=6.8036989E+17
  
  !do n = 1,num_scalar(ngrid)
  !srf - somente para gases, veja rotina original para todos os demais escalares
  do n=num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED)+1,num_scalar(ngrid)
     !------------------- TMP (SO para 1 primeiro gas)
     if(n>num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED)+1) CYCLE ! TMP
     !------------------- TMP
  
  
     call calc_masscon(mzp,mxp,myp,scalar_tab(n,ngrid)%var_p &
  		 ,grid_g(ngrid)%dxt    (1,1)	&
  		 ,grid_g(ngrid)%dyt    (1,1)	&
  		 ,basic_g(ngrid)%dn0 &
  		 , masscon)
  
  print*,'N masscon=',n, masscon,100.*(masscon-initalmass)/initalmass

  enddo
  end SUBROUTINE check_mass_conservation
  !-------------------------------------------------------

  SUBROUTINE calc_masscon(m1,m2,m3,var_p, dxt,dyt,dn0,masscon)
  use mem_grid, only : dzt
  implicit none
  integer, intent(in) :: m1,m2,m3
  real, dimension(m1,m2,m3), intent(in) :: var_p,dn0
  real, intent(out) :: masscon
  real, dimension(m2,m3), intent(in) :: dxt,dyt
  integer :: i,j,k
  real vol
  
  masscon=0.
  do j=2,m3-1; do i=2,m2-1
  
  
      do k=2,m1-1
    !print*,'dzt =', 1./dzt(k),1./dxt(i,j),1./dyt(i,j)
    vol= 1./( dzt(k)*dxt(i,j)*dyt(i,j)) 
    masscon=masscon+var_p(k,i,j)*dn0(k,i,j)*vol

  enddo;enddo;enddo

  end SUBROUTINE calc_masscon
  !-------------------------------------------------------
  !-------------------------------------------------------
END MODULE monotonic_adv
